Divergent responses of butterflies and bees to burning and grazing management in tallgrass prairies

Abstract Butterflies and bees contribute significantly to grassland biodiversity and play important roles as pollinators and herbivores. Grassland conservation and management must be seen through the lens of insect conservation and management if these species are to thrive. In North America, grasslands are a product of climate and natural disturbances such as fire and grazing. These natural disturbances have changed considerably since European colonization and subsequent landscape fragmentation. The aim of this study was to better understand the impacts of fire and grazing management on butterfly and bee communities in tallgrass prairie, enabling land managers and conservationists to better protect and manage remnant prairie. We examined butterfly and bee abundance, species richness, and diversity in Minnesota tallgrass prairies managed by grazing or fire. In 2016 and 2017, we surveyed butterflies, bees, vegetation, and surrounding land use at 20 remnant prairies (10 burned and 10 grazed) with known management histories. Butterfly and bee abundance at our study sites were significantly negatively correlated. Butterfly abundance, but not species richness, was higher in burned than grazed prairies, and prairie‐associated grass‐feeding butterflies were more abundant at sites with higher plant species richness. Bee abundance was unrelated to management type but was higher at sites with sandier soils; bee species richness was positively associated with forb frequency. These findings highlight the challenges of designing management plans tailored to wide groups of pollinators and the potential pitfalls of using one group of pollinators as indicators for another. They also point to the importance of a mosaic of management practices across the prairie landscape.


| INTRODUC TI ON
Thoughtful and informed land management is necessary if declining native grasslands and their inhabitants are to persist. Butterflies and bees contribute significantly to grassland biodiversity and play important roles in ecosystem functioning. Butterfly adults are incidental pollinators, butterfly larvae are important herbivores (Scoble, 1992), and all life stages serve as food sources for birds and other animals. Bees are considered the most important pollinators both globally and in tallgrass prairie (Grimaldi, 1999). Worldwide declines in insect diversity and abundance are increasingly welldocumented (e.g., Biesmeijer et al., 2006;Cameron et al., 2011;Wagner, 2020), including the butterflies and bees that are the subjects of this study. Prairie specialist butterflies are declining in Iowa, Wisconsin, Minnesota, and Illinois tallgrass prairies (e.g., Schlicht et al., 2009;Swengel & Swengel, 2013;Swengel et al., 2011), and 10 of the 15 endangered, threatened, or special concern butterfly species in Minnesota are associated with tallgrass prairie (Minnesota Department of Natural Resources, 2013). Some of these, including the threatened Dakota Skipper (Hesperia dacotae) and federally endangered Poweshiek Skipperling (Oarisma poweshiek), were once among the most common butterflies in tallgrass prairie (Dana, 1991;Schlicht & Orwig, 1998). The federally endangered rusty-patched bumble bee (Bombus affinis), which occurs in Minnesota, was abundant only twenty years ago and is now rarely found across most of its historic range (USFWS, 2016). The imperiled status of these insects warns us that common species are not resistant to declines faced by insects as a whole. Extensive changes to natural disturbance regimes in the Minnesota tallgrass prairie, coupled with habitat loss and fragmentation, are potential drivers of declines of once ubiquitous insect species. It is therefore increasingly important that grassland conservation and management take insects into consideration when developing management plans.
North American prairie evolved and was maintained for tens of thousands of years through ungulate grazing, lightning-ignited fires, and indigenous fire management (Anderson, 2006;Middleton, 2013), which reduced woody plant growth. Land managers often attempt to mimic natural fire and grazing disturbances through prescribed fire and cattle grazing management (Brudvig et al., 2007). However, with so much of the historic extent of prairie gone and what remains scattered across a fragmented landscape, managers face increasing challenges when seeking to maintain remnant prairie (prairie that has never been plowed or converted to agriculture). At least 98% of Minnesota's approx. 7,285,000 hectares of the tallgrass prairie has been converted to agriculture or otherwise lost, and other tallgrass prairie states have suffered similar losses (Samson et al., 2004). This habitat loss and fragmentation results in substantial threats to biodiversity (e.g., Brudvig et al., 2015;Fahrig, 2003;Haddad et al., 2015;Summerville & Crist, 2001).
Although fire and grazing occurred concurrently or in response to one another historically (Anderson, 2006), managers today are often faced with the choice of either burning or grazing based on logistic (e.g., having the infrastructure to manage cattle or sufficient distance from human habitation to apply fire) or economic (e.g., willing livestock owners to graze on the remnant prairie or available trained personnel to apply fire) feasibility. Fire management has become more challenging as prairie remnants become fragmented, smaller, and more isolated. Managers are often constrained by the increased presence of humans, farmland, and roadways in the landscape because they must account for wind direction and smoke and the risk of fire escaping (Toledo et al., 2013). Additionally, leaving unburned refugia for prairie obligate insects (Swengel et al., 2011) becomes more difficult in smaller remnants. Although spatially dependent, these constraints can result in fire frequencies that are lower than many resource managers would consider optimal, and also lower than are used in most research studies on fire effects (e.g., Collins & Calabrese, 2012;Dickson et al., 2019). Management must respond to local conditions, and Midwestern tallgrass prairies rarely, if ever, receive the frequent fire that is more typical in places like Konza Prairie, where much of the influential research on fire and grazing originated.
Conservation grazing, in which domestic herbivores are used to further conservation goals (Asensio & Lauenroth, 2012), is one way to reduce potential threats of fire. However, today's conservation grazing is done almost exclusively with domesticated cattle, which preferentially graze different vegetation, prefer wetter areas, and move with different herd patterns than bison (Allred et al., 2011;Kohl et al., 2013;Plumb & Dodd, 1993). Grazing also requires partnerships with livestock owners who support conservation outcomes, and the additional fencing and water infrastructure required often makes grazing impractical. In addition to logistical challenges, it is not always clear which management strategy will produce the desired ecological outcomes. Grasslands are disturbance-dependent landscapes, but there remains much debate about how best to practice disturbance management in the current landscape, especially with regard to insect conservation (e.g., Buckles & Harmon-Threatt, 2019;Henderson et al., 2018).
Studies examining the impacts of fire and grazing management on butterflies and bees often find inconsistent results. Panzer (2002) and Thom et al. (2015) report that prairie remnant-dependent butterfly species that overwinter above ground as eggs, larvae, or pupae are particularly vulnerable to fire, especially if there are few nearby refugia from which butterflies may recolonize a burned site (Driscoll et al., 2010;Swengel & Swengel, 2007). Swengel (1998) found that, in general, the majority of butterfly species studied occurred in greater abundance under mowing and grazing management than under rotational-burning management. On the other hand, butterflies typically absent during the time when fires are set, such as monarchs (Danaus plexippus) (Leone et al., 2019;Moranz et al., 2012) and other migratory species, or that are in life stages that occur underground (e.g., Maculinea spp in Europe (Nowicki et al., 2015)) may not suffer negative effects of burning but instead benefit from habitat improvement. Vogel et al. (2007) found that while butterfly species richness did not differ between management practices, butterfly diversity indices were highest in burn-only sites and species composition differed by management. In comparison, bees' responses to fire or grazing are influenced by their life history, including nesting location. Those that nest 10 cm or deeper underground (75% of ground-nesting taxa) tolerate most grassland fires, which typically do not raise soil temperatures to lethal levels nor for lethal durations (Cane & Neff, 2011;DeBano et al., 1998). Fires can be more dangerous for insects that nest aboveground due to both nest combustion and lethal temperatures (Tooker & Hanks, 2004). Results have been mixed regarding the impact of grazing on grassland bees. Kimoto et al. (2012) found that grazing intensity had no significant effect on total bee abundance or species richness in the central Oregon prairie. There were differences in response between genera, with greater intensity grazing more negatively impacting Bombus (bumble bee) than Lasioglossum (sweat bee) abundance. Increased grazing intensity was also associated with a lower Shannon diversity in bees in the early season, potentially due to declines in floral resources.
However, Carvell (2001) found a greater abundance of bumble bees in pastures grazed by cattle within the past year.
Butterflies are sometimes used as pollinator "indicator" taxa in ecological studies (Thomas, 2005), due to the comparative ease of sampling and identifying butterflies compared with bees.
However, there is debate about their usefulness as indicators. Davis et al. (2008) found that butterfly and bee diversity were negatively correlated in Iowa tallgrass prairies, although management practices were not considered in their study. Management plans that assume similar responses from different pollinator groups may only benefit some species, while others are left out. It is essential for grassland management and butterfly and bee conservation that these assumptions are tested.
To inform better management of tallgrass prairie butterflies and bees, we investigated how bees and butterflies respond after ≥11 years of fire or grazing management as practiced by resource managers. We thus are considering the cumulative effects over time of these management practices on bees and butterflies, rather than the direct and immediate effects of fire or grazing on the organisms.
Our goals were to assess (1) the effects of conservation grazing versus prescribed fire management on butterfly and bee abundance and richness and (2) whether butterflies and bees differ in their responses to fire versus grazing management. Specifically, we investigated the abundance and species richness of all observed butterflies and bees, as well as subsets of each: resource-user butterflies, which represent observed butterflies seen using resources within managed sites, as opposed to flying through; prairie-associated grass-feeding butterflies, which we were interested in because of their relation to species of conservation concern; and soil-excavating groundnesting bees, which are among the most abundant and speciose bee taxa typically collected in bee bowls.
While site management is important in shaping prairie bee and butterfly communities, it does not occur in isolation. We hypothesized that both butterfly and bee communities would be affected by management practices, but that their responses to fire vs grazing would differ and be mediated by local and landscape factors such as patch size, prairie habitat availability in the landscape, floral and host plant resources, and soil texture. Habitat patch size and the amount of suitable habitat in the surrounding landscape are known to positively impact bee and butterfly communities (Denning & Foster, 2018;Robinson et al., 2014;Topp et al., 2021). We expected these to be positively associated with the abundance or diversity of butterflies and bees. Forbs provide nectar for bees and butterflies, and pollen for bees (Denning & Foster, 2018;Öckinger & Smith, 2006;Winfree et al., 2011). We thus expected the abundance of butterflies and bees to be positively associated with forb frequency. Host plant resource availability is important in shaping butterfly communities (Dennis et al., 2011). Nine of Minnesota's endangered, threatened, and at-risk butterfly species feed on native graminoids, as do all other members of the subfamily Hesperiinae (Hesperiidae) (Narem & Meyer, 2017;Scott, 1986). We expected butterfly and bee diversity to be positively associated with plant species richness and prairie-associated grass-feeding butterfly diversity to be positively associated with native graminoid frequency.
For ground-nesting bees, soil accessibility (i.e., bare soil) and texture are vitally important. Fire initially increases bare soil exposure, which can provide ground-nesting bees with more nesting opportunities.
Grazing also increases the amount of bare soil in grasslands, with higher intensities resulting in more bare ground (Kimoto et al., 2012).
Despite the increase in bare soil exposure, we also expected groundnesting bee abundance to be negatively associated with grazing frequency and intensity, as cows can lead to increased soil compaction and inundation of soils (Alaoui et al., 2018;Batey, 2009;Buckles & Harmon-Threatt, 2019).
Because the effects of land management can take years to appear, and because we wanted to provide insights directly relevant to the types of prairies with which land managers work in Minnesota, we chose study sites that were managed at least once during the eleven years prior to this study by state, federal, and private land managers and that were exclusively burned or grazed for at least eleven years prior to the beginning of our study.

| Study sites
We chose 20 remnant, tallgrass prairie sites within the prairie parkland province in Minnesota ( Figure 1) from candidate sites that had all been exclusively either burned or grazed by cattle between 2005 and 2015 (10 burned, 10 grazed). Sites represented a range of sizes (1.13-144.7 ha), prairie habitat in the surrounding landscape (0.15%-68%), years managed (1-13 years), time since fire (2-9 years), and cattle stocking rates (0.17-2.9 AUM, Animal Unit Month) (Appendix 1, Table A1). Management records, permits, and permissions were granted by owners (the US Fish and Wildlife Service, Minnesota Department of Natural Resources (MN DNR), The Nature Conservancy, and private landowners).
We created a 1.5-km buffer (Greenleaf et al., 2007;Lane et al., 2020) around each site and calculated the percent of the prairie surrounding each site, not including the site itself, using ArcMap (v 10.5.1); see full methods in Larson et al. (2018).
Briefly, we calculated the percent of land classified as prairie within the CropScape cropland data layer (Han et al., 2014), MN DNR native prairie and Reinvest in Minnesota-MN Geospatial data (MN DNR), and South Dakota State University potentially undisturbed land (Bauman et al., 2016) within the 1.5-km buffer around each site.

| Sampling methods
Butterflies and bees were surveyed at sites three times in both 2016 (June 15 to August 31) and 2017 (May 14 to August 18), for a total of 117 surveys. One grazed site, (G-1) was only surveyed in 2017.
To address phenology differences across the north-south range of sites, we surveyed sites from south to north. We conducted 60 bee and butterfly surveys at burned sites and 57 at grazed sites during the study.
To minimize the effect of time of day on sampling, sites surveyed in the afternoon during one visit were surveyed in the morning during the next visit and vice versa. To reduce weatherrelated sampling variability, insect surveys were conducted between 09:30 h and 18:30 h (with 2 exceptions when surveying finished between 18:30 h and 19:00 h), when temperatures were above 20°C, sustained winds were less than 20 km/h, and cloud cover was <70% (15 exceptions) with no precipitation (Moranz et al., 2014;Pollard & Yates, 1993). Using available soil drainage data, we delineated wet, mesic, and/or dry prairie polygons for each site. Transects were delineated within each prairie-type polygon prior to field sampling, and oriented parallel to elevation gradients. The total length of insect transects was the same at all sites: Butterfly transects were 400-m long and bee transects were 180-m long, sharing the same beginning points. For insect transect survey purposes, we sampled each prairie type in proportion to its portion of the total site area. For example, if a site was delineated as 50% mesic prairie, 40% wet prairie, and 10% dry prairie, we conducted 200 m of butterfly transects surveys along transects in the mesic prairie, 160 m in wet prairie, and 40 m in dry prairie. If a site was 20% mesic prairie and 80% wet prairie, we would conduct 80 m of butterfly transect surveys in mesic prairie and 320 m in wet prairie. We similarly distributed bee bowls proportionally along wet, mesic, and/or dry prairie transects. At some sites, one continuous transect did not fit and transects were broken into smaller sections due to prairie type, shape, or size; at F I G U R E 1 Map of burned (B 1-10; triangles) and grazed (G 1-10; circles) tallgrass prairie study sites within the prairie parkland province in Minnesota.
Potential Tallgrass Prairie Extent these sites, transects were at least 20-m apart to avoid sampling redundancy.

| Butterfly surveys and identification
All butterfly surveys were conducted by the same observer using two methods. First, we used a modification of the standardized Pollard Walk for relative abundance (e.g., Pollard, 1977;Pollard & Yates, 1993), during which we walked transects at a steady pace of 10 m/minute and recorded each individual butterfly seen within a 5-m imaginary box in front of the observer: 2.5 m on each side, 5 m ahead, and up to 5 m above the ground. This method provides relative abundance data and is used in the analyses that follow.

| Bee surveys and identification
During each site visit, bees were surveyed in two ways, passively via pan traps ("bee bowls"), and actively, via netting, to achieve the most complete account of species at the sites. All site visits also included a time-constrained meandering walk in which bees were netted when observed on flowers; the timer was stopped while processing bees. The length of the meandering walk is scaled with site size, lasting between 30 min and 2 h. Netted bees were placed individually in a glassine envelope, labeled with the date, time, and site name, and kept in an ethyl acetate jar until frozen for later processing. Data from this method were only used to supplement species richness data from bee bowls but not used in abundance analyses. All bee identification took place in the laboratory with the use of a stereomicroscope, using the following keys and guides: Gibbs (2010Gibbs ( , 2011, Gibbs et al. (2013), LaBerge (1967, 1980, Laverty and Harder (1988), Ribble (1968), and Williams et al. (2014).
Discover Life (Ascher & Pickering, 2018) was also consulted. A table of all bee species identified is included in Appendix 2, Table A3. (invasive graminoid) frequencies were calculated based on the presence of each detected species (number of occupied plots/total plots) (Appendix 1). We used plant frequency because sampling occurred throughout the growing season, so cover in early surveys would not be comparable to cover in later surveys (Elzinga et al., 1998). Poa pratensis and Bromus inermis are invasive thatch-forming grasses

| Vegetation surveys
that land managers seek to control through fire and grazing management. Five 10-cm × 2.54-cm soil cores were collected at each site along a randomly selected vegetation transect once in either 2016 or 2017, from which the proportion of sand was calculated (Appendix 1). Vegetation and soil methods are described fully in Larson et al. (2020).

| Butterfly response variables
Four measures of butterfly abundance were modeled separately: total butterfly abundance, resource-user abundance, non-monarch butterfly abundance, and prairie-associated grass-feeding butterfly abundance. We summed total butterfly abundances from three survey visits at each site in 2016 and 2017 separately for each year, resulting in an index of butterfly relative abundance, which we hereafter refer to as total butterfly abundance (n = 39). We analyzed a subset of total abundance that included only butterflies that were observed using resources at a site (i.e., we removed butterflies from the analyses that had only been observed flying and not nectaring, basking, mating, ovipositing, or performing other activities related to site resource-use), hereafter, "resource-users." We did this because butterflies only observed flying at a site may not be impacted by local management, especially for smaller sites and larger, more mobile species. We analyzed a subset of total abundance that included only non-monarch butterflies. We did this because monarchs accounted for a large proportion of butterfly observations at our sites (Table A2) and we previously found them to be positively associated with fire at the same sites (Leone et al., 2019). We also analyzed prairie-associated butterflies whose larvae feed on grasses based on Schlicht et al. (2007), Narem and Meyer (2017), and personal communications with local butterfly experts. These species, listed in Table A2, are of interest to us because many prairie-associated butterflies in Minnesota and the Midwestern United States feed exclusively on native grasses in their larval stages, including the once-common but now federally endangered Oarisma poweshiek and federally threatened Hesperia dacotae. The abundances of many species within the prairie-associated group were too low to allow for species-specific analyses, so we grouped all prairie-associated grassfeeding species together for analyses.
Butterfly species richness was estimated using the Chao 2 estimator (Chao, 1984;Colwell & Coddington, 1994). The suite of species observed at a site can be sensitive to bias due to the size of the site, the conditions during site visits, and the effort during surveys (Chao et al., 2014). Observed species richness can thus be an unreliable measure of the full community at a site, especially considering that some species are very rare and therefore unlikely to be detected. We calculated Chao 2 as: The term S Y,T represents the number of species observed during transect surveys plus meandering walk surveys at site T during year Y, L Y,T is the number of species that occur in only one sample from site T during year Y, and M Y,T is the number of species that occur in exactly two samples at site T during year Y. The estimated richness and the observed richness become more similar as the ratio of unique species to doubly observed species gets smaller. This is based upon the assumption that in the true community, many fewer species should occur in a single sample than in two samples. Thus, as the ratio of L to M gets smaller, the Chao 2 estimator approaches S. As species richness is a count of discrete species, a Poisson distribution is appropriate for models. We rounded the Chao 2 estimator to the nearest integer and used the fossil package (Vavrek, 2011) in R 3.6.2 (R Core Team, 2019) to calculate this estimator for each site in 2016 and 2017. Hereafter, "butterfly species richness" refers to the Chao 2 estimated value.

| Bee response variables
Bee abundance was adjusted to account for the loss of bee bowls at grazed sites when cattle were present. The adjusted bee abundance was calculated as: This calculation estimates the number of bees that would have been collected had an entire set of traps (30 bowls × 3 visits = 90) been recovered in a given year. Rounding to the nearest integer allows for the use of the Poisson distribution, which is appropriate for count data. For most site visits, where all 30 bee bowls were recovered, the adjusted bee abundance and raw bee abundance were identical. Hereafter, "bee abundance" will refer to adjusted bee abundance.
Bee species richness was estimated using the Chao 2 estimator described above. Hereafter, "bee species richness" refers to the Chao 2 estimated value. We also analyzed a subset of total abundance and bee species richness that included only bees that excavate nests underground. Ground-nesting bees were categorized according to an in-progress database from Bartomeus et al. (2013).

| Butterfly and bee models
The response variables described above were analyzed using Poisson distributed generalized linear mixed-effects models (GLMMs).
Predictor variables were selected a priori based on the literature and included management type as a categorical variable (burned, grazed), the percent of prairie within 1.5 km, site area, forb frequency, and the combined frequency of two invasive, thatch-forming graminoids (Poa pratensis and Bromus inermis). Butterfly models included plant species richness and native graminoid frequency, to account for potential host plant associations. Bee models included the proportion within site, were included as random effects in all models. We tested the likelihood ratio between models with the random effects structure of year nested within site vs. models with just site as a random effect. We found that models that included year within site differed significantly from models that included only site as a random effect, indicating that these models can parameterize temporal variation despite the grouping factor having only two levels (Gomes, 2022).
This method accounts for the well-documented phenomenon of interannual variation in insect pollinators (e.g., Herrera, 1988;Price et al., 2005). We report the random intercept variance values for the final models in Appendix 3, Table A4.
We did not include additional management variables in our models because they were associated with management type (stocking rate was only relevant at grazed sites, time since fire only relevant at burned sites, and number of years managed not comparable between burned and grazed sites [Appendix 1, Table A1]). Instead, we built GLMMs with subsets of the data (burn-only sites and graze-only sites) to examine associations between all response variables and the predictor variables stocking rate and number of years managed at grazed sites and time since fire and number of years managed at burned sites.
We compared adjusted abundance and species richness for butterflies and bees using the Spearman's rank correlation and the function cor.test from the stats package in R (R Core Team, 2019).
For both butterflies and bees, analyses were conducted in R 3.6.2 (R Core Team, 2019) using the glmer function from the lme4 package (Bates et al., 2015)  Total butterfly abundance was close to two times higher at sites managed with fire than those managed with grazing (z = −2.332, p = .0197; Figure 2a); all other predictor variables were removed during backward elimination. The abundance of butterfly resourceusers and non-monarch butterflies was also higher at burned sites than grazed sites (z = −2.22, p = .0264; Figure 2b and z = −0.4177, p = .0413, respectively), with management type as the only significant predictor variable after backward elimination in both cases. The abundance of prairie-associated grass-feeding butterflies was similar in burned and grazed sites (z = 0.069, p = .9448). The model with the lowest AIC value (ΔAIC > 2; Arnold, 2010) for prairie-associated grass-feeding butterfly abundance after backward elimination included only plant species richness (z = 1.680, p = .0929), which was positively, but not significantly, associated with abundance (alpha = 0.05). Other habitat variables did not explain any variation in the number of prairie-associated grass-feeding butterflies. We observed no prairie-associated grass-feeding butterflies in either 2016 or 2017 in four sites, two additional sites had no observations in 2016, and fewer than five individuals were observed at four of the occupied sites.
F I G U R E 2 Total (a) and resource-user (b) butterfly abundance at burned (B) and grazed (G) tallgrass prairie sites within the prairie parkland province in Minnesota, USA. Box plots depict the minimum, first quartile, median, third quartile, and maximum, with outliers depicted as single points.
Total, resource-user, non-monarch, and prairie-associated grassfeeding butterfly abundance was similar at grazed sites with different stocking rates and numbers of years managed and at burned sites with different times since fire and number of years managed.
Univariate model results for all butterfly abundance response variables are presented in Appendix 4, Tables A5-A7.

| Butterfly species richness
We observed 39 butterfly species over the course of two summers; 36 in 2016 and 32 in 2017; 34 at sites managed with fire and 34 at sites managed with grazing (

| Butterfly community composition
The first two axes in the butterfly NMS (stress = 9.5 with 45 itera- Speyeria cybele (r = .584) were most strongly negatively associated with years managed (Appendix 5, Table A13).

| Bee abundance
We collected 11,969 bees in bowl traps in the summers of 2016 and 2017. A univariate analysis of the effect of the duration of bee bowl deployment on adjusted bee abundance showed no significant cor- Total bee abundance was higher at sites with sandier soils (z = 2.421, p = .0155; Figure 4a); no other variables were significant in the final multivariate model. The abundance of soil-excavating ground-nesting bees was also higher at sites with sandier soils (z = 2.456, p = .014).
Neither time since fire nor the number of years managed with fire had a significant effect on total bee abundance or soil-excavating F I G U R E 3 Sites within-species space for nonmetric multidimensional scaling analysis of butterflies on grazed (pink squares) and burned (blue triangles) sites, axes 1 and 2. The second axis represented 28% of the variation in the data and was correlated with years managed (r = −.511).
Vectors are proportional to the strength of the correlation with the axes. See Appendix 5, Table A13 for all correlations between butterfly species and NMS axes.
ground-nesting bee abundance, nor did stocking rate or the number of years managed with grazing. Univariate model results for bee abundance response variables are presented in Appendix 4, Tables A9-A12.

| Bee species richness
We identified 119 species (30 genera) in our 2016 and 2017 collections. Sixty-two specimens were not identified as species or species complex and were not included in richness analyses. One hundred two species were collected at burned sites, 25 of which were exclusive to burned sites, and 94 species at grazed sites, 17 of which were exclusive to grazed sites (Table A3). Twenty-seven species were represented by only a single specimen ("singletons"), and 18 species were represented by two specimens ("doubletons") ( are soil-excavating ground-nesters and 11 (9.2%) occupy existing cavities (Table A3). Approximately 88% of individuals collected (11,004 of 12,540) are in the family Halictidae, bees that are mostly small ground-nesters that generally prefer sandier soils (Cane, 1991;Potts & Willmer, 1997).
The final multivariate model for bee species richness included forb frequency, which was positively associated with species richness (z = 2.99, p = .0028; Figure 4b), and site area, which was negatively, but not significantly, correlated with species richness (z = −1.511, p = .1308).
None of the predictor variables tested were associated with ground-nesting bee species richness.
Neither time since fire nor the number of years managed with fire had a significant effect on total bee species richness or soilexcavating ground-nesting bee species richness. The number of years managed with grazing had a significant effect on total bee species richness (z = −2.367, p = .018), with fewer bee species found at sites grazed more frequently. Neither the stocking rate nor the number of years managed with grazing had a significant effect on soil-excavating ground-nesting bee species richness.
Univariate model results for bee species richness are presented in Appendix 4.

| Bee community composition
The  Table A14).

| DISCUSS ION
Butterfly abundance differed between burned and grazed remnant prairie, but bee abundance and species richness were related to sand and forb frequency at our study sites. These findings highlight the F I G U R E 4 Relationship between (a) bee abundance and proportion sand and (b) bee species richness (Chao2) and forb frequency at sites in 2016 (black) and 2017 (white).
challenges of designing coherent management plans tailored to wide groups of pollinators and the dangers of using one group of pollinators as indicators for another ( Table 1).
We expected any butterfly response to management to be me- Our finding that butterfly species richness did not differ based on management is consistent with others (Moranz et al., 2012;Vogel et al., 2007). However, our finding of higher butterfly abundance at sites managed with fire compared with grazing is more nuanced, and previous studies are more varied in their results. We note that fire frequencies at our sites (1-3 times in a 11-year period) were sometimes much lower than in otherwise comparable studies. Vogel et al. (2007) found that most habitat generalists did not differ in abundance among management practices, although they reported a higher abundance of D. plexippus and Colias eurytheme in sites managed with only grazing compared to those managed with only burning; burn frequencies varied from 1-3 times in an 8-year period.
In comparison, we found similar C. eurytheme abundance between grazed and burned sites and nearly three times as many D. plexippus at burned compared with grazed sites (Table A2), which may be driving some of the patterns in overall abundance in our models (see Leone et al. (2019) for a more in-depth analysis of D. plexippus). Our finding that Speyeria idalia abundance was higher at burned compared with grazed sites is supported by Vogel et al. (2007). Our findings are also consistent with Moranz et al. (2012), who reported the highest population densities of C. pegala, S. idalia, and D. plexippus in burn-only treatments (Table A2). By contrast, Vogel et al., 2007 found that among habitat specialists, Cercyonis pegala abundance was higher in grazed than burned sites. Clearly, species identities influence butterfly responses to management.
In contrast to our results for total and resource-user butterfly abundance, the abundance of prairie-associated grass-feeding butterflies did not differ between burned and grazed sites in our study.
The positive relationship between the abundance of these butterflies and plant species richness combined with the fact that plant species richness did not differ between management types at our study sites (Larson et al., 2020) suggests that this association is unrelated to fire or grazing. Many grass-feeding prairie-associated butterfly species have seen precipitous declines in recent decades; in fact, many such species were not observed during this study  Schlicht et al., 2009;Swengel et al., 2011). The species in this group that we did observe were generally in low abundances. However, community composition and NMS results help differentiate species responses. Of the five species we included in the prairieassociated grass-feeding butterfly group, only C. pegala, the most abundant species in this group, was more abundant at burned sites than grazed sites ( Table A2). Three of the remaining four species, H. F I G U R E 5 Sites within species space for nonmetric multidimensional scaling analysis of bees on grazed (pink squares) and burned (blue triangles) sites, axes 1 and 2. The first axis represented 52% of the variation in the data and was correlated with proportion sand (r = −.722) and proportion clay (r = .514).
Vectors are proportional to the strength of the correlation with the axes. See Appendix 5, Table A14 for all correlations between bee species and NMS axes.
leonardus, P. themistocles, and C. tullia were observed only at grazed sites (Table A2) and had strong positive associations in NMS with years managed (Figure 3; Appendix 5). Although many of these butterflies are included in studies of tallgrass prairie butterflies (e.g., Davis et al., 2008;Moranz et al., 2012;Schlicht et al., 2009;Swengel et al., 2011;Vogel et al., 2007), few studies have compared the impacts of management strategies for them (but see Swengel, 1998).
Low abundances, specialized life histories, and association with plant species richness suggest that a more targeted study may be needed for these species of concern.
Differences in time since fire have been found to influence butterfly abundance. However, we found no effect of time since fire on butterfly abundance or on butterfly species richness in our study.
This is in contrast to Vogel et al. (2010), who reported a positive association between butterfly abundance and time since fire, with 50-to 70-month recovery times postfire in the Loess Hills of Iowa.
Significant positive postburn responses to fire have also been documented for monarch butterflies and their milkweed host plants within one to two years following fire (e.g., Baum & Sharber, 2012;Rudolph et al., 2006). By contrast, lower butterfly abundance has been documented at burn-only prairies than burn-and-graze prairies with a fire rotation of 2-6 years (Vogel et al., 2007). Because none of our sites were burned during the study or the preceding year (2015), differences in butterfly abundance are unlikely to reflect qualitative differences in nectar or host plant resources as a direct result of fire.
Butterfly populations could have recovered from any negative impacts of fire at our study sites prior to surveying.
Another possible explanation for the higher abundance of butterflies at burned prairies compared with grazed prairies is a negative effect of grazing, rather than a positive association with fire.

TA B L E 1
Butterfly and bee responses to fire versus grazing in Minnesota tallgrass prairie: significant associations between response variables and predictor variables in final models after backward selection. Note: Positive association with management type indicates that values were higher in burned (vs. grazed) sites and "ns" is not significant. Significance codes: 0.001 "***"; 0.01 "**"; 0.05 "*". number of years a site was grazed were not correlated with butterfly abundance in this study, there may be indirect effects of grazing that we did not quantify; we only measured the frequency and not the quality of plant species. Although forb frequency did not differ between our burned and grazed study sites (Larson et al., 2020), frequent fire has been shown to increase nectar availability (Rudolph et al., 2006); grazing may also reduce the amount of floral resources.

Management type
We did not quantify floral resources but did observe cattle consuming flowers. Grazing reduces vegetation height, and several studies have found that butterflies prefer taller vegetation (Berg et al., 2013;Öckinger & Smith, 2006;Poyry et al., 2006).
Neither bee abundance nor species richness were influenced by management type in our study; this does not necessarily mean that bees do not respond to management but may mean that burning and Comparisons to historical collections would be a worthy avenue for future research. The few bee species restricted to burned or grazed sites (Table A3) are only represented by one or two individuals, making any conclusions about their true exclusivity impossible. These species may be rarely occurring, or rarely captured using our techniques, making their detection at either management type just as unlikely.
Bees generally and the subset of soil-excavating ground-nesters were more abundant in sites with sandier soils. Different bee species prefer to nest in soils of different textures, although relatively few bees are associated with clay-rich soils; most prefer sandy loams (Cane, 1991). These soils are easier to excavate and less susceptible to flooding than silt-or clay-heavy soils (Skiba & Ball, 2002). Our bee community analyses also support the importance of soil texture in shaping the bee community; the proportions of sand and clay in soils were relatively strongly correlated with the first axis of the NMS, which explained most of the variation in the community ( Figure 5).
Analysis of soil-excavating ground-nesting bees, which represent the most abundantly collected bees in our samples, showed the sandiness of soils as the only significant predictor of their abundance. This indicates that soil-excavating ground-nesters are driving patterns of bee abundance. It may also indicate, as noted below, that biases in taxa collected by bee bowls are influencing analyses.
The response of the bee community to grazing is not a simple one. While soil-excavating ground-nesters have an important influence on models of total bee abundance, there are also signals of the importance of aboveground nesters and nonexcavators in grazed prairies. The frequency of grazing, measured as the number of years within the previous 10 years that a prairie was grazed, had a significant negative effect on total bee species richness. Like Kimoto et al. (2012), our best fit model did not include stocking rate as a significant predictor of bee abundance or species richness. Kimoto et al. (2012), offers us another point of comparison; they found that the abundance of the generally soil-excavating genus Lasioglossum was less negatively impacted by grazing than the generally nonexcavating genus Bombus. Contrary to our expectations, there was a negative relationship between grazing frequency and species richness of the whole prairie bee community while the community of soil-excavating ground-nesters was not impacted in our study. While we expected that increased grazing frequency would compact soils, making bee nests more prone to inundation (Alaoui et al., 2018;Batey, 2009) and thus limiting soil-excavating ground-nesting bees' ability to make use of grazed sites (Buckles & Harmon-Threatt, 2019), we detected nothing to indicate this.
Although the frequency of grazing had a significant effect on total bee species richness at grazed sites, forb frequency was the only significant predictor of bee species richness across all sites, a finding in line with prior research that documented floral resource availability as a limiting factor for bees (e.g., Inari et al., 2012;Ogilvie et al., 2017;Roulston & Goodell, 2011) and other pollinator communities (Sjödin, 2007). At our study sites, forb frequency itself was not significantly impacted by management type (Larson et al., 2020), but the lack of an association between management type and bee species richness was surprising, nonetheless.  (2018) was not explained by a shifting floral community; rather the same plant species seen at unburned sites flowered longer at burned sites.
We found that an 11-year history of burning and grazing, in isolation, does not predict bee abundance or species richness.
This equivalency of abundance and richness between burned and grazed prairies, as well as the lack of significant distinction between the communities making up burned and grazed prairies in which semi-natural grassland habitats exist has no significant effect on bee abundance. Importantly, the majority of individuals we collected do not rely solely upon prairie fragments; the four most abundant bee species in our study, Lasioglossum pruinosum, L. albipenne, L. versatum, and Augochlorella aurata (54% of all bees collected) are widely distributed across North America in various habitats (Coelho, 2004;Gibbs, 2011). Management effects, while potentially destructive to some individuals, may have minimal impacts on these species at a landscape scale, leading to seeming equivalency between management approaches for bee abundance and species richness. This is supported by our finding that management did not affect the abundance or presence of soil-excavating ground-nesters, which include the four species listed above.
Butterfly and bee abundance were negatively correlated, and we found no correlation between butterfly and bee species richness.
Because associated predictor variables differed between butterflies and bees in our models, we urge caution in the use of one as an indicator of habitat suitability for the other. Davis et al. (2008) (Thomas, 2005), our study highlights their inadequacy as predictors of bee abundance and richness.
This retrospective study offered a duration of a single management type that would have been impossible to achieve through experimental manipulation in the time frame of this project.
Additionally, the tallgrass prairie is a rare resource, and land managers are tasked with protecting and promoting that resource, whether for the public or for their herds. An observational study allowed us to work with land managers without compromising their missions; many of the managers worked with us in the hope that our findings could inform future management decisions on these very same lands. However, the retrospective nature of this study imposes some limitations. The lack of experimental manipulation made parsing out the direct and indirect effects of management difficult. We were also limited in our ability to control the extent of variation in factors unrelated to management, such as site area or latitude. Controls for variation in sites had to be made at the time of site selection, but it is possible that variables outside of our consideration, such as site history before 2005, could obscure signals. Interannual variation in insects, including among bees and butterflies, is well-documented (e.g., Fishbein & Venable, 1996;Herrera, 1988;Price et al., 2005).
We did our best to account for this background temporal variation by including year as a random effect in our models and reporting the random intercept variance in Appendix 3. However, we recognize the limitation imposed by two years of sampling data across highly variable populations.
While bee bowls have been widely used in recent decades, they were not used historically, making comparisons with previous indices of prairie bee communities difficult .
Bee bowl samples are biased towards certain taxonomic groups, with members of the family Halictidae over-represented as compared to other collection means (e.g., Droege et al., 2010;Geroff et al., 2014;Griffin et al., 2017). Bees may also be drawn to bee bowls from the surrounding areas, especially when flowers are scarce (Kuhlman et al., 2021), making our samples a measure of both the surveyed site and the surrounding matrix of grassland, agriculture, and development (Baum & Wallen, 2011;Roulston et al., 2007). These last two points-the taxonomic bias and the potential attraction outside of the study area-may be driving results. The effect of sandy soils may be amplified by the fact that bee bowls attract the very bees that prefer sandy soils, obscuring other signals. The lack of significance of management type may be because bee bowl samples are drawing bees in from the wider area, where disturbance and habitats are more homogenized.
However, we did attempt to curb these limitations by including a meandering walk to capture a wider breadth of bee species richness than found in bee bowls alone. Additionally, our analyses included the percent of prairie in the 1.5 km surrounding the sites, thus providing a measure of the broader habitat matrix that could account for unknown variation brought by bee bowls' attraction of bees from outside the study site. Ultimately, collection methods will always shape the sample of the community they provide. We present these limitations here in acknowledgment of that fact and encourage future studies to take them into account.
Patch-burn grazing, in which cattle are set to graze on recently burned vegetation, is increasingly implemented to create a patchwork of heavily and lightly disturbed areas (Fuhlendorf et al., 2009;Helzer & Steuter, 2005) and can thus promote diverse and heterogeneous plant communities. The extent to which this creates good bee or butterfly habitat is unclear, however (

| CON CLUS IONS
The fact that bee and butterfly communities, with the exception of butterfly abundance, did not differ between sites managed with grazing or infrequent fire over a 13-year period can be taken as an encouraging sign; the management practice that is most appropriate and practical in a given situation can be used without concern about harming invertebrate communities broadly, although some species appear to do better under one management practice. Burning, at least at sites managed with fire 1-3 times over 11 years, appears to support higher butterfly abundance, although this may be the result of untested variables and not the direct result of fire. Some species are more likely to be found in grazed sites and species composition differs with the number of years a site is managed. A variety of management strategies across sites is therefore important to support the entire suite of bee and butterfly species. Managers interested in promoting bee abundance and diversity might consider increasing forb frequency and targeting sites with sandier soils for acquisition, preservation, or future restoration. data curation (supporting); formal analysis (supporting); funding acquisition (lead); investigation (equal); methodology (equal); project administration (equal); resources (equal); supervision (lead); validation (supporting); visualization (supporting); writing -original draft (supporting); writing -review and editing (equal).

ACK N OWLED G M ENTS
We thank the Minnesota Environment and Natural Resources

CO N FLI C T O F I NTE R E S T
All authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data and metadata will be accessible through the Data Repository   , 2005-2015. "Mgmt" stands for Management, "Pct" stands for Percent (0%-100% prairie in surrounding 1.5 km site buffer), and "Prop" stands for Proportion (0-1 proportion sand within each site). Vegetation and soil variables were rounded to two decimal places.

TA B L E A 1 (Continued)
TA B L E A 3 Adjusted abundance and presence of bee species observed in 2016 and 2017 at burned and grazed tallgrass prairie sites within the prairie parkland province in Minnesota, USA. Note: ΔAIC indicates the relative differences between each univariate model and the "best-ranked" (minAIC) model.

TA B L E A 7
Univariate model results showing the response of prairie-associated grass-feeding butterfly abundance to each predictor variable tested. Note: ΔAIC indicates the relative differences between each univariate model and the "best-ranked" (minAIC) model.

A PPE N D I X 5
Pearson and Kendall Correlations with ordination axes from nonmetric multidimensional scaling analysis of butterfly and bee species.

TA B L E A 1 3 Pearson and Kendall
Correlations with ordination axes from nonmetric multidimensional scaling analysis of butterfly species.